2-2

(({{(Subscript[a, 1]), (Subscript[a, 2])}, {(Subscript[b, 1]), (Subscript[b, 2])}}))

Construction & application of a Random Transition Matrix

The elements of such a matrix are "transition probabilities." Writing

𝕋=(({{(Subscript[T, 1,1]), (Subscript[T, 1,2]), (Subscript[T, 1,3]), (Subscript[T, 1,4])}, {(Subscript[T, 2,1]), (Subscript[T, 2,2]), (Subscript[T, 2,3]), (Subscript[T, 2,4])}, {(Subscript[T, 3,1]), (Subscript[T, 3,2]), (Subscript[T, 3,3]), (Subscript[T, 3,4])}, {(Subscript[T, 4,1]), (Subscript[T, 4,2]), (Subscript[T, 4,3]), (Subscript[T, 4,4])}}));

we understand

(Subscript[T, j,k])=probability of transition j←k  from site k to site j

(({{(Subscript[p, 1])[n+1]}, {(Subscript[p, 2])[n+1]}, {(Subscript[p, 3])[n+1]}, {(Subscript[p, 4])[n+1]}}))=𝕋.(({{(Subscript[p, 1])[n]}, {(Subscript[p, 2])[n]}, {(Subscript[p, 3])[n]}, {(Subscript[p, 4])[n]}}))

where

(Subscript[p, k])[n]= probability that the "walker"is at site k after n steps

Necessarily, the elements on any row must sum to unity:

(∑_(j=1)^N)(Subscript[p, j,k])=1

To construct a random instance of such a matrix (in the case N = 5) I construct a square array of numbers all of which (in this instance) fall randomly on the interval [0, 1]:

RowData=Table[Table[RandomReal[],{j,1,5}],{k,1,5}]

Normalizing so that the elements on each row sum to unity, I construct the random transition matrix

𝕋=Transpose[Table[((RowData[[k]])/(Total[RowData[[k]]])),{k,1,5}]]

and check that the elements of each column do in fact sum to unity:

Table[Total[Transpose[𝕋][[k]]],{k,1,5}]

𝕋//MatrixForm

Random walks are deterministic processes, so we expect them to be reversible. Introducing

𝕌=Inverse[𝕋];

𝕌//MatrixForm

we see that many of the elements of the inverted transition matrix cannot be interpreted as "transition probabilities," (some fall outside the interval [0, 1], some even wear negative signs) though they do—for non-obvious reasons—conform to one of the characteristic features of transition matrices: the elements of each column sum to unity:

Table[Total[Transpose[𝕌][[k]]],{k,1,5}]

Looking to the eigenvalues of 𝕋

λ[k_]:=Eigenvalues[𝕋][[k]]

Table[{k,λ[k]},{k,1,5}]//TableForm

which Mathematica presents in order of descending norm

Table[Norm[λ[k]],{k,1,5}]

of which one is very much larger than the others. We observe that some eigenvalues are real, some complex, and that the complex eigenvalues occur in conjugate pairs. Those reality/complixity/conjugacy distinctions are repeated when one looks to the corresponding eigenvectors:

v[k_]:=Transpose[{Eigenvectors[𝕋][[k]]}]

Table[{k,MatrixForm[v[k]]},{k,1,5}]//TableForm

We check that we do in fact have now in hand 5 solutions of the "eigen-equation:"

Table[Chop[𝕋.v[k]-λ[k]v[k]]//MatrixForm,{k,1,5}]

The eigenvectors are, as supplied by Mathematica, normalized

Table[Norm[v[k]],{k,1,5}]

but are not orthogonal (we had no reason to expect them to be):

Table[Chop[Conjugate[Transpose[v[j]]].v[k]],{j,1,5},{k,1,5}]//MatrixForm

It is a curious fact that all eigenvectors except the leading eigenvector have elements that sum to zero:

Table[Chop[Total[Transpose[v[k]][[1]]]],{k,1,5}]

The leading eigenvector is unique in that it can be "unitized," meaning multiplied by a factor so as to cause the elements to sum to unity:

V=((v[1])/(Total[Transpose[v[1]][[1]]]));

V//MatrixForm
Total[Transpose[V][[1]]]

Notice that V possesses the structure of a probability vector.

𝕋.(({{1}, {0}, {0}, {0}, {0}}))//MatrixForm

Generally,

𝕋.(probability vector)=another probability vector

which I illustrate by simple example:

𝕋.(({{1}, {0}, {0}, {0}, {0}}))//MatrixForm
Total[Transpose[%][[1]]]

Iteration, whatever the initial probability vector, leads inevitably back to V:

?MatrixPower

MatrixForm[MatrixPower[𝕋,10].(({{1}, {0}, {0}, {0}, {0}}))-V]
Total[Transpose[%][[1]]]

MatrixForm[MatrixPower[𝕋,20].(({{1}, {0}, {0}, {0}, {0}}))-V]
Total[Transpose[%][[1]]]

To illustrate the point more generally, I construct a random initial probability vector

x=Table[RandomReal[],{k,1,5}]

(Subscript[p, 0])=Transpose[{(x/(Total[x]))}];
(Subscript[p, 0])//MatrixForm

and obtain

MatrixForm[MatrixPower[𝕋,10].(Subscript[p, 0])-V]
Total[Transpose[%][[1]]]

The process

(Subscript[p, 0])→→→⋯→V

is deterministic, and reversible-in-principle

MatrixPower[𝕋,-10].MatrixPower[𝕋,10].(Subscript[p, 0])-(Subscript[p, 0])//MatrixForm

but if the number of iterative steps is too large the reversal is sabotaged by numerical imprecision:

MatrixPower[𝕋,-20].MatrixPower[𝕋,20].(Subscript[p, 0])-(Subscript[p, 0])//MatrixForm

Total[Transpose[%][[1]]]

Repetition with a higher-dimensional example

BigRowData = Table[Table[RandomReal[], {j, 1, 10}], {k, 1, 10}];

𝕋𝕋=Transpose[Table[((BigRowData[[k]])/(Total[BigRowData[[k]]])),{k,1,10}]];

Table[Total[Transpose[𝕋𝕋][[k]]],{k,1,10}]

λλ[k_]:=Eigenvalues[𝕋𝕋][[k]]

Table[{k,λλ[k]},{k,1,10}]//TableForm

Table[Norm[λλ[k]],{k,1,10}]

ListPlot[%]

Notice that again, one of the norms is very much larger than  the others.

vv[k_]:=Transpose[{Eigenvectors[𝕋𝕋][[k]]}]

Table[Chop[𝕋𝕋.vv[k]-λλ[k]vv[k]]//MatrixForm,{k,1,10}]

Table[Norm[vv[k]],{k,1,10}]

Table[Chop[Total[Transpose[vv[k]][[1]]]],{k,1,10}]

VV=((vv[1])/(Total[Transpose[vv[1]][[1]]]));
VV//MatrixForm
Total[Transpose[VV][[1]]]

xx=Table[RandomReal[],{k,1,10}];
(Subscript[pp, 0])=Transpose[{(xx/(Total[xx]))}];
(Subscript[pp, 0])//MatrixForm
Total[Transpose[(Subscript[pp, 0])][[1]]]

I demonstrate convergence to VV as the number of iterations increases

Table[{2k,MatrixForm[MatrixPower[𝕋𝕋,2k].(Subscript[pp, 0])-VV]},{k,1,8}]//TableForm

and (working at default machine precision) the computational breakdown of reversibility after about 8 iterations:

Table[{2k,MatrixPower[𝕋𝕋,-2k].MatrixPower[𝕋𝕋,2k].(Subscript[pp, 0])-(Subscript[pp, 0])//MatrixForm},{k,1,8}]//TableForm

Imposition of detailed balancing

The transition matrices considered thus far are asymmetric:

(Subscript[p, j,k])≠(Subscript[p, k,j])

We now impose the principle that the probability of jumping j→k is the same as the probability of jumping j←k :

(Subscript[p, j,k])=(Subscript[p, k,j])

which renders 𝕋 symmetric.

Construction technique:

The method I will use to construct symmetric transition matrices is sketched below:

(({{(Subscript[A, 1]), (Subscript[A, 2]), (Subscript[A, 3]), (Subscript[A, 4]), (Subscript[A, 5])}, {□, □, □, □, □}, {□, □, □, □, □}, {□, □, □, □, □}, {□, □, □, □, □}})): arrange to have ∑(Subscript[A, k])=1
(({{(Subscript[A, 1]), (Subscript[A, 2]), (Subscript[A, 3]), (Subscript[A, 4]), (Subscript[A, 5])}, {(Subscript[A, 2]), (Subscript[B, 1]), (Subscript[B, 2]), (Subscript[B, 3]), (Subscript[B, 4])}, {(Subscript[A, 3]), □, □, □, □}, {(Subscript[A, 4]), □, □, □, □}, {(Subscript[A, 5]), □, □, □, □}})): arrange to have (Subscript[A, 2])+∑(Subscript[B, k])=1
(({{(Subscript[A, 1]), (Subscript[A, 2]), (Subscript[A, 3]), (Subscript[A, 4]), (Subscript[A, 5])}, {(Subscript[A, 2]), (Subscript[B, 1]), (Subscript[B, 2]), (Subscript[B, 3]), (Subscript[B, 4])}, {(Subscript[A, 3]), (Subscript[B, 2]), (Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3])}, {(Subscript[A, 4]), (Subscript[B, 3]), □, □, □}, {(Subscript[A, 5]), (Subscript[B, 4]), □, □, □}})): arrange to have (Subscript[A, 3])+(Subscript[B, 2])+∑(Subscript[C, k])=1
   etc.

I implement the procedure (which will be found to work except when it fails to work, due to the intrusion of minus signs, from which it provides no protection), proceeding on the assumption that matrix elements are generated randomly by flat distribution on the unit interval, and use the implementation to construct an 8-element population of such 5×5 symmetric transition matrices:

a=Table[RandomReal[],{k,1,5}];
b=Table[RandomReal[],{k,1,4}];
c=Table[RandomReal[],{k,1,3}];
d=Table[RandomReal[],{k,1,2}];
e=Table[RandomReal[],{k,1,1}];
(Subscript[Row, 1])=(a/(Total[a]));
(Subscript[Row, 2])=Join[{(Subscript[Row, 1])[[2]]},((1-(Subscript[Row, 1])[[2]])/(Total[b]))b];
(Subscript[Row, 3])=Join[{(Subscript[Row, 1])[[3]],(Subscript[Row, 2])[[3]]},((1-(Subscript[Row, 1])[[3]]-(Subscript[Row, 2])[[3]])/(Total[c]))c];
(Subscript[Row, 4])=Join[{(Subscript[Row, 1])[[4]],(Subscript[Row, 2])[[4]],(Subscript[Row, 3])[[4]]},((1-(Subscript[Row, 1])[[4]]-(Subscript[Row, 2])[[4]]-(Subscript[Row, 3])[[4]])/(Total[d]))d];
(Subscript[Row, 5])=Join[{(Subscript[Row, 1])[[5]],(Subscript[Row, 2])[[5]],(Subscript[Row, 3])[[5]],(Subscript[Row, 4])[[5]]},((1-(Subscript[Row, 1])[[5]]-(Subscript[Row, 2])[[5]]-(Subscript[Row, 3])[[5]]-(Subscript[Row, 4])[[5]])/(Total[e]))e];
(Subscript[𝕊, 8])={(Subscript[Row, 1]),(Subscript[Row, 2]),(Subscript[Row, 3]),(Subscript[Row, 4]),(Subscript[Row, 5])};

Table[{k,MatrixForm[(Subscript[𝕊, k])]},{k,1,8}]

We verify that (as is evident to the eye) each of those matrices is symmetric

Table[(Subscript[𝕊, k])==Transpose[(Subscript[𝕊, k])],{k,1,8}]

and that the rows/columns of each sum to unity:

TableForm[Table[{j,{Table[Total[(Subscript[𝕊, j])[[k]]],{k,1,5}]}},{j,1,8}]]

Notice, however, that (in this particular sample) (Subscript[𝕊, 5]), (Subscript[𝕊, 6]) and (Subscript[𝕊, 8]) must be discarded because they contain elements that wear minus signs, and therefore do not admit of interpretation as transition probabilities.

Eigenvalues

From real symmetry it follows familiarly that all eigenvalues are real.

EigenvalueTable=Table[Eigenvalues[(Subscript[𝕊, k])],{k,1,8}]
ListPlot[Flatten[EigenvalueTable]]

Interestingly, the spectral value 1 is universal,  and all other eigenvalues have absolute value < 1:

ListPlot[Abs[Flatten[EigenvalueTable]]]

Eigenvectors

Looking to the respective populations of eigenvectors

Flatten[Table[{j,TableForm[Table[MatrixForm[Transpose[{Eigenvectors[(Subscript[𝕊, j])][[k]]}]],{k,1,5}]]},{j,1,8}]]

we note that the leading eigenvector (the one associated with the leading eigenvalue λ = 1) is the same in every instance, that all of its elements are equal, and that all wear the same sign (which could without loss of generality be taken to be—i.e., adjusted to be—positive). 

The real symmetry of the 𝕊-matrices causes their eigenvectors (which Mathematica delivers in normalized form) to be orthogonal, which we verify:

Table[{n,Table[Chop[{Eigenvectors[(Subscript[𝕊, n])][[j]]}.Transpose[{Eigenvectors[(Subscript[𝕊, n])][[k]]}]][[1]][[1]],{j,1,5},{k,1,5}]//TableForm},{n,1,8}]

From the fact that the 𝕊-matrices possess the structure peculiar to transition matrices it appears to follow that the elements of all but the leading eigenvectors sum to zero:

Table[Chop[Total[Eigenvectors[(Subscript[𝕊, j])][[k]]]],{j,1,8},{k,1,5}]

The leading eigenvectors—and only the leading eigenvectors—can be scaled to become "probability vectors:"

V[j_]:=Transpose[{((Eigenvectors[(Subscript[𝕊, j])][[1]])/(Total[Eigenvectors[(Subscript[𝕊, j])][[1]]]))}];

Table[MatrixForm[V[j]],{j,1,8}]

Those vectors—the same in every instance—refer to a flat or "equal likelyhood" distribution, an "equipartition of probability." Interestingly, these results hold even when occasional minus signs intrude into the structure of 𝕊.

Interation and its asymptotic endstate

(Subscript[𝕊, 1]).(({{1}, {0}, {0}, {0}, {0}}))//MatrixForm

Table[MatrixForm[MatrixPower[(Subscript[𝕊, 1]),n].(({{1}, {0}, {0}, {0}, {0}}))],{n,0,10}]

p=Table[RandomReal[],{k,1,5}]

(Subscript[p, 0])=Transpose[{(p/(Total[p]))}];
(Subscript[p, 0])//MatrixForm

Table[MatrixForm[MatrixPower[(Subscript[𝕊, 1]),n].(Subscript[p, 0])],{n,0,10}]

P[n_]:=Transpose[MatrixPower[(Subscript[𝕊, 1]),n].(Subscript[p, 0])][[1]]

MovingElement[k_]:=Table[P[n][[k]],{n,0,10}]

ListLinePlot[{MovingElement[1],MovingElement[2],MovingElement[3],MovingElement[4],MovingElement[5]},PlotRange->{.1,.35}]

Table[Chop[MatrixForm[MatrixPower[(Subscript[𝕊, 1]),-n].MatrixPower[(Subscript[𝕊, 1]),n]]],{n,0,15}]

MatrixPower[(Subscript[𝕊, 1]),12]//MatrixForm

Det[%]//Chop

Inverse[(Subscript[𝕊, 1])]//MatrixForm

Inverse[(Subscript[𝕊, 1])]

Table[Total[Inverse[(Subscript[𝕊, 1])][[k]]],{k,1,5}]